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We present a study of dynamic properties of liquid aluminum using density- functional theory within the 
local-density (LDA) and generalized gradient (GGA) approximations. We determine the temperature 
dependence of the self-diffusion coefficient as well the viscosity using direct methods. Comparisons with 
experimental data favor the LDA approximation to compute dynamic properties of liquid aluminum. We 
show that the GGA approximation induce more important backscattering effects due to an enhancement of 
the icosahedral short range order (ISRO) that impact directly dynamic properties like the self- diffusion 
coefficient. All these results are then used to test the Stokes-Einstein relation and the universal scaling law 
relating the diffusion coefficient and the excess entropy of a liquid. 

I ncreasing demands to optimize methods and facilities for material processing technologies such as melting, 
I refining, or casting metallic alloys are boosting the precise knowledge of thermophysical properties of these 
I alloys. The main goals are an improvement of the final product quality, an enhancement of the process 
efficiency, and an economical consumption of resources and energy. More particularly, as the solidification 
process of a liquid alloy has a strong impact on the structure and properties of the solid material, the under- 
standing of the behavior and more particularly the knowledge of the thermophysical properties of molten alloys 
prior to solidification are essential for the development of materials with predetermined characteristics. For 
instance, it is known that properties of lightweight aluminum-based alloys are affected by the temperature 
conditions of melting and casting 1 , and the knowledge of the diffusion process in their liquid phases is hence a 
necessity. 

Among all thermophysical properties, diffusion coefficients are known to play a prominent role for a detailed 
understanding of the solidification process including crystal growth 2 " 4 and vitrification 5 . Quite surprisingly, 
experimental values for liquid aluminum and its alloys are scarce. Their direct measurements in liquid alloys 
are mostly made using so-called capillary tube methods that require isotopes as tracers. However, the influences of 
convective flow on the diffusion profile during annealing can be a big problem 6 . As a consequence, self- diffusion 
coefficients D can be overestimated and published values often differ by a factor of 2. Moreover, radioactive 
isotopes are not available for aluminum. D are also often obtained indirectly from experimental determinations of 
the shear viscosity coefficients, r\, easier to measure and the use of the Stokes-Einstein (SE) relation 7 , D SE = k B T/ 
2n Rrj, k B being Boltzmann's constant and R is an effective atomic diameter taken as the position of the first peak 
of the pair-correlation function. If the applicability of the SE relation to liquid metals is well accepted above their 
melting temperatures, the scatter of experimental viscosity values is important in the case of aluminum 8 , reflecting 
difficulty of obtaining accurate data. This can be due to the fact that Al has a high affinity to oxygen. Solid oxide 
particles can possibly be trapped between the melt and the crucible, affecting drastically the flow profile and 
therefore the viscosity. For instance, the scatter between measurements made with different techniques is ±14% 
in the case of aluminum while it is only ±6% for iron 8 . 

Recently, it has been shown for some liquid metals like Ni, Ti, or Cu 9-11 that self- diffusion coefficients can be 
obtained accurately from the use of quasielastic neutron scattering (QNS) in a large temperature range, via the 
determination of the incoherent intermediate scattering function. For liquid Al, the situation is more complicated 
since its incoherent cross section is two orders of magnitude lower than that of liquid Cu. Very recently, Demmel 
et al. 12 and Kargl et al. 13 report temperature- dependent self- diffusion coefficients of liquid Al, with different 
assumptions to extract the incoherent contribution from the scattering function. The two sets of experimental 
data are in good agreement and display an Arrhenius-type behavior. The absolute values of D at and around the 
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melting temperature agree also well with those obtained from the 
Stokes-Einstein relation. However, at high temperatures, i.e. 1200 K, 
the discrepancy becomes important 12 . 

From a theoretical point of view, molecular dynamics (MD) is 
a widely used simulation technique for investigating the structural 
and dynamic properties in the liquid state 14 . In classical molecular 
dynamics, the simulated particles are moved using forces extracted 
from semi- empirical interatomic potentials, electronic degrees of 
freedom remaining hidden. By contrast, ab initio molecular dynamics 
(AIMD) within the density functional theory (DFT) are based on 
quantum-many-body forces computed via the Hellmann-Feynman 
theorem, which is straightforward using plane-wave basis 15 . 

Properties of liquid Al have been studied extensively using the semi- 
empirical embedded-atom method (EAM) formalism. However Becker 
et al 16 have shown that diffusion coefficients are strongly dependent on 
the different implementations of the EAM potential. Moreover, the 
best agreement leads to an underestimation of the experimental values 
that becomes more important at high temperatures. 

Within the DFT framework, Alfe and Gillan 17 have shown that the 
calculated self- diffusion and shear viscosity coefficients are sensitive 
to the choice of the density. Their calculations gave D values within 
the range 0.52-0.68 A 2 /ps while n values are in the range 1.4- 
2.2 mPa.s. Gonzalez et al. 18 used an orbital-free DFT method to 
compute dynamic properties of liquid Al at T = 943 and 1323 K. 
This approach underestimates the experimental data, the discrep- 
ancy becoming very important at high temperatures. Such a result 
shows that the accuracy of the orbital-free DFT approach to compute 
dynamic properties of liquid Al is likely to be limited by the use of an 
approximate electronic kinetic energy functional and a local pseu- 
dopotential. In a more recent work, the same group 19 performed 
Kohn-Sham ab initio molecular dynamics simulations and reported 
the same diffusion coefficient D as in Ref. 17. They also computed the 
shear viscosity from the transverse current correlation function. 
They obtained n = 1.05 mPa.s, which is smaller than the value 
obtained by Alfe and Gillan 17 , using the definition of rj as the time 
integral of the auto- correlation function of the off-diagonal compo- 
nents of the stress tensor. 

In this paper, we address the important issue of the determination 
of the dynamic properties of liquid Al using ab initio molecular 
dynamics simulations. Both the local- density approximation 20 ' 21 
(LDA) and the generalized gradient approximation (GGA) in the 
Perdew-Burke-Ernzerhof 22 (PBE) form for the exchange- correlation 
energy were investigated. We show that the determination of the self- 
diffusion coefficient as well as the viscosity depends on the exchange- 
correlation functional used in DFT calculations. More particularly, 
GGA based calculations of the diffusion coefficient underestimates 
both LDA based values and experimental data. We demonstrate that 
this behavior is related to the detailed description of the short-range 
order and more particularly to the enhancement of the icosahedral 
short-range order using the GGA approximation. All these results 
are then used to test the Stokes-Einstein relation and the universal 
scaling law relating the diffusion coefficient and the excess entropy of 
a liquid 23 ' 24 . 

Results 

In this section, we determine the single atom as well collective 
dynamics that give access to transport properties such as the self- 
diffusion coefficient, D, and viscosity, rj. 

Self-diffusion coefficient. In a first step, we consider the single-atom 
dynamics through the velocity auto -correlation function, which is 
defined as: 

^(0=(v,(0-v,(f = 0)) (1) 

where v,-(f) are the position and velocity of atom i at time U and the 
angular brackets denote the average over time origins as well as 
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Figure 1 | Evolution of the self-diffusion coefficients of liquid Al as a 
function of inverse temperature. The solid circles are AIMD results with 
LDA and the open ones are those with GGA. The dashed lines are their 
respective Arrhenius fit. The open triangle and the star correspond to 
previous ab initio calculations with LDA 9,11 . The open diamonds and full 
squares are experimental data of Refs. 6 and 7, respectively, along with their 
Arrhenius fit (dash-dotted lines). 

atoms. The time integration of the velocity auto -correlation 
function gives access to the self- diffusion coefficient, D. This is 
equivalent to the determination of D using the linear slope at long 
times of the mean-square displacement 7,25 , and require shorter 
simulation times to have statistically meaningful results. 

The temperature dependence of D using the LDA and GGA 
approximations is reported in Fig. 1. We can see that GGA calcula- 
tions underestimate LDA calculations by about 20% above the melt- 
ing point and by about 30% in the undercooled state. LDA 
calculations are closer to the experimental values 1213 and similar to 
previous LDA- AIMD simulations at T = 1000 K 1719 . Assuming an 
Arrhenius -type behavior for the diffusion process in the liquid state, 
an activation energy E can be derived. As shown in Fig. 1, the fit 
describes the data quite well in the whole range of temperature. The 
derived activation energy is E = 249 ± 5 and 260 ± 5 meV respect- 
ively for the LDA and GGA approximations. Activation energies 
have been also obtained from experimental data. The reported values 
are 280 ± 70 (Ref. 13) and 274 ± 4 meV (Ref. 12) and are only 
slightly larger than our calculated values. 

Viscosity. Another interesting dynamical magnitude is the 
transverse current correlation function C T (q,t) which has the 
advantage of yielding a generalized g- dependent shear viscosity 
from which the hydrodynamic limit can be evaluated 26 ' 27 . It is a 
straightforward approach based on the decay of the correlation of 
the atom motions, which is directly linked to the viscosity, and does 
not require the use of the forces as needed in the Green-Kubo 
formalism 28 . We briefly describe here the main points of the 
formalism. Firstly, the function Cj(q,f) is defined as 

C T (q,t)=^(r T (q,t)J T (q,t)), (2) 

where Ji{q,t) represents the transverse current that is expressed along 
the x direction as 

N 

h(q,t) = V| >W ex Pfe(0] • (3) 

i = 1 

The quantity v iyX (t) is the ^-component of the velocity of atom i and q 
is a wave vector along the z direction, recalling that q = (2n/L)(n x , n y , 
n z ) are reciprocal vectors which have to be compatible with the length 
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L = V 1/3 of the simulation cell, and n x , n y and n z are integers. Two 
formally identical expressions are also used for y and z directions. 
The zero limit Laplace transform Or(g,£ = 0) gives rise to the in- 
dependent shear viscosity 

pk B T 



mq 2 CT(q,z = 0) 



(4) 



k B and m stand respectively for Boltzmann's constant and mass of the 
atoms. In the hydrodynamic limit, i.e. by extrapolating rj(q) to q = 0, 
the shear viscosity rj can be obtained. The extrapolation is performed 
with the function proposed by Alley and Alder 26 

^=TTW (5) 

which is used in general to represent the viscosity of a dense hard- 
sphere packing and has proven to work well for liquid Ni 29 . Fig. 2 



shows that the fitting gives a good representation of the LDA and 
GGA results of g- dependent viscosity for liquid Al at different 
temperatures. 

At 1000 K, the LDA calculated value is 1.17 mPa.s close to the 
assessed value 2 of 1.18 mPa.s. The values given by previous AIMD 
simulations using LDA are 1.4 and 1.05 mPa.s, respectively obtained 
from the time integral of the auto-correlation function of the off- 
diagonal components of the stress tensor 17 and from the transverse 
current correlation function 19 . The value afforded by the Stokes- 
Einstein relation is 1.21 mPa.s, by using R = 2.72 A, which corre- 
sponds to the first peak of the LDA pair-correlation function. The 
good agreement between this value and the value of 1.17 mPa.s 
calculated above supports the validity of using the SE relation to 
describe dynamic properties of liquid Al using LDA. For GGA, the 
value obtained from the SE relation is 1.53 mPa.s (R = 2.72 A cor- 
responding to the first peak of the GGA pair-correlation function) 
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Figure 2 | g-dependent viscosity of liquid Al obtained from our AIMD with (a) LDA, and (b) GGA. The symbols are the calculations using the transverse 
current correlation function for T= 875 K or 850 K (open diamonds), T= 1000 K (open squares), T= 1125 K (open triangles), and T= 1250 K (open 
circles). The solid lines are the fits using Eq. (5). 
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Figure 3 | Viscosity of liquid Al as a function of temperature. The full circles and full squares correspond to our AIMD results, respectively with 
LDA and GGA, as obtained with the transverse current correlation functions. The open circles and open squares correspond to viscosities inferred from 
the self- diffusion coefficients, respectively with LDA and GGA, determined from our AIMD results using the Stokes-Einstein relation. The open 
triangle and the star correspond to previous ab initio calculations 911 . The dashed, dash-dotted, and dash-dot-dotted lines correspond to the experimental 
values of Ref. 23 and 24, and the assessment of Ref. 2, respectively. 



which can be compared to that obtained by the direct method, 
namely 1.37 mPa.s. Both values overestimate the assessed value 
and their difference is larger than that obtained in LDA. This also 
holds for the other temperatures, as displayed in Fig. 3. It can be seen 
that the results of LDA are close to the assessed values of Ref. 8. 
Moreover, we note that they are also close to the experimental values 
of Iida and Shirashi 30 and remain consistent with those of Battezzati 
and Greer 31 , both experimental data sets being retained in the com- 
pilation by Dinsdale and Quested 32 . 



Discussion 

Differences between GGA and LDA calculations of the self- diffusion 
coefficient and the viscosity are significant, more than 15% for the 
viscosity and 20% for diffusion. To understand these discrepancies, 
we inspect some computed structural properties using both GGA 
and LDA functionals. 

In Fig. 4, we display the pair distribution function g(r) of liquid Al 
for all the thermodynamic states as described in the preceding sec- 
tion. We can see that our simulations agree well with experimental 




r(A) 

Figure 4 | Evolution of the pair-correlation function with temperature. The open circles and triangles correspond to experimental values of Refs. 26 
(1023 K) and 27 (1123 K and 1223 K), respectively. The curves for 1125 K, 1000 K and 875 K are shifted upwards by an amount of 1,2 and 3, respectively. 
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Figure 5 | (a) Evolution of the coordination number with temperature. The solid squares are the results of LDA and the open ones are those of GGA. The 
open triangles are the corrected experimental values of Ref. 18; (b) Evolution of abundances of the most important pairs for liquid Al with temperature. 
The solid symbols are the results of LDA and the open ones are those of GGA. Error bars are typically of the order of 0.01. 



X-ray diffraction measurements 33,34 as well as with previous ab initio 
simulations 19 for the liquid state. We do not see any appreciable 
difference between the two sets of calculations but it is known that 
changes which may occur in the short ranger order (SRO) require 
more detailed analysis 35,36 . 

Integration of the computed pair distribution function g(r) up to 
its first minimum gives the coordination number, which is displayed 
in Fig. 5(a) as function of temperature. We find an increase of coor- 
dination number from N c = 11.3 ±0.1 at T= 1250 KtoiV c = 11.8 ± 
0.1 at T = 875 K while GGA calculations give N c = 11.6 ± 0.1 at T = 
1250 K and N c = 12.1 ± 0.1 at 850 K. We obtain that LDA calcula- 
tions give coordination numbers smaller than those obtained using 
GGA calculations but differences remain small and cannot explain 



differences obtained in the dynamic properties. Note that the experi- 
mental values of N c reported by Mauro et al 34 increase from 12.9 ± 
0.2 to 13.1 ± 0.2 from T = 1123 K to T = 1273 K. However the 
authors of Ref. 33 estimated that these values are overestimated by 
about 8%. Correction of the experimental coordination numbers by 
this amount would lead to values in good agreement with our AIMD 
results and with a similar temperature evolution as shown in 
Fig. 5(a). 

More insight into the structural changes is gained by analyzing the 
inherent structures and characterizing the local environment sur- 
rounding each atomic pair that contributes to the first peak of g(r) 
in terms of the number and properties of common nearest neighbors 
of the pair under consideration 35 . In Fig. 5(b), we report the most 
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Figure 6 | Velocity auto-correlation function obtained for liquid Al at 1000 K with LDA (solid lines) and GGA (dashed lines). Inset: same comparison 
in the vicinity of the main minimum for 875 or 850 K, 1000 K, 1 125 K, and 1250 K. For the last three temperatures the curves are shifted upwards by an 
amount of 0.05, 0.1 and 0.15, respectively. 



abundant pairs, averaged over the ten inherent configurations for 
each temperature and for both LDA and GGA calculation, i.e., 142 X 
(sum of 1422 and 1421), 1431 and 15XX (sum of 1551 and 1541) 
pairs found in both calculations. The number of 15 X X bonded pairs 
is a direct measure of the degree of icosahedral ordering including 
both perfect and distorted icosahedral motifs while 142X bonded 
pairs are characteristic of close packed structures. The 1431 pairs 
either can be considered as distorted icosahedra or distorted close- 
packed structures 37 " 40 . As they are similar in both approximations 
and hold steady with temperature, they are not considered to be 
responsible for differences between dynamic properties using either 
GGA or LDA. The same comment holds for the 142 X pairs. On the 
contrary, Fig. 5(b) shows that the number of 15 XX pairs is more 
important using the GGA approximation. Our findings indicate that 
the icosahedral short range order is more pronounced when using 
the GGA approximation. We can notice also that the icosahedral 
15 XX fraction always exceeds the close-packed 142 X fraction and 
it is the only fraction that displays a linear increase as the temperature 
decreases. This trend is in fine agreement with results based on the 
atomistic cluster alignment method 41 . 

To correlate ISRO with dynamic properties, we focus on the back- 
scattering regime that is predominant for pure liquid Al because of its 
high density. The backflow induced by a moving atom increases the 
probability of an atom to jump back toward its initial position. It is 
characterized by a well-pronounced minimum in the velocity auto- 
correlation function which is very sensitive to the local liquid density. 
In Fig. 6, we report the velocity auto-correlation function computed 
at T = 1000 K using the LDA and GGA approximations. The first 
minimum becomes deeper using GGA, which reveals an increase of 
backscattering effects. Note that this difference becomes more 
important when the temperature decreases, as shown in the inset. 
This variation can be related to the increase of the icosahedral motifs 
using the GGA approximation, as discussed in Ref. 42. Icosahedral 
motifs are known to be the most compact local structures and con- 
sequently the backscattering is more pronounced in the liquid which 
presents the highest ISRO and the diffusivity is smaller. 

Understanding the connection between the transport coefficients 
and structural properties of dense liquids can be also tackled using 



the correspondence between the transport coefficients and the 
internal entropy established by Rosenfeld 24 or Dzugutov 23 . The latter 
author proposed a universal scaling relationship between the excess 
entropy of a liquid, S E , and a dimensionless form of the diffusion 
coefficient, D*. The excess entropy is the total entropy of the liquid 
minus the ideal gas contribution. The relationship can be written as: 



Where D* is given by: 



D*= 0.049 exp(S £ ) 



D*=DT~ l a~ 



(6) 



(?) 



In Eq. (7), D is the diffusion coefficient and T is the Enskog 
collision frequency which can be calculated for the temperature T 
and the density p as: 



T = 4a 2 g((T)p 



nk B T\ 1/2 



(8) 
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Figure 7 | The scaled diffusion coefficient vs the S 2 approximation using 
ab initio calculations and vs the excess entropy calculated by the 
Carnahan Starling (CS) equation. 
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Table I | values of the packing fraction £, and isothermal compres- 
sibility xt as a function of temperature T 



T(K) Z = npa 3 /6 



Xt~- 



(4 



1 - Z) 4 /{2pk B K 
{) + (!- fl 4 ) 
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where m and a are the atomic mass and the hard-sphere diameter. 
g(<j) represents the value of the pair correlation function g(r) at the 
contact distance. 

In addition, Dzugutov assumed that the excess entropy can be 
approximated by the two -body expression denoted: 



2np\ (g(r)\n\g(r)}-\g(r)-l})r 2 dr 



(9) 



If the S 2 description is a very good approximation for simple model 
systems like Lennard- Jones systems, Hoyt et al 43 have shown that it 
is less accurate for the density dependent many-body EAM poten- 
tials. Therefore, we test Dzugutov's idea using our computed 
diffusivities and radial distribution functions within the LDA 
approximation. Results displayed in Fig. 7 suggest that the S 2 
approximation is also not completely reliable when utilizing either 
LDA or GGA calculations. 

Thus, as for embedded atom potentials, a more reliable determina- 
tion of the excess entropy is needed. The excess entropy can be 
expressed formally as a multiparticle correlation expansion: S E = 
S 2 + S 3 + . . ., where S n is the entropy contribution due to n -particle 
spatial correlations 44 . Using this approach is not obvious since 
Yokoyama et al 45 have shown that the expansion converges very 
slowly for metallic systems since the four-body and even higher order 
terms of correlations have non-negligible contributions. Hoyt et al 43 
used a lambda integration technique to get the excess entropy at T = 
1500 K followed by Monte Carlo simulations to evaluate the excess 
entropy at any desired temperature. 

Here we adopt a different strategy by making use of the Carnahan - 
Starling equation which is known to give the best analytical descrip- 
tion of the thermodynamics of hard-sphere systems as long as the 
packing density cj is smaller than 0.5 (Ref. 46). Note that the packing 
density is given by £ = npo 3 l6. In this case, the excess entropy can be 
derived from the Carnahan- Starling equation as 



= g(3g-4) 2(2-0 JH_ 
1 (1-0 2 (1-0 3 W 



(10) 



All entropy calculations can be done using the packing fraction £, and 
its evolution with temperature. We determine £, using the experi- 
mental densities 8 and a values obtained by considering that the pair- 
correlation entropy of the hard sphere model, S^ s , is equal to the 
pair- correlation entropy provided by LDA calculations, S\ DA . The 
pair correlation function of the hard sphere model, ^ s (r), is com- 
puted by the thermodynamically self- consistent integral equation 
method 47 . The evolution of the packing fraction with temperature 
is given only by the evolution of the experimental density with tem- 
perature since ^ DA {r) and then a do not vary with temperature as 
obtained in our simulations. We can mention that cj varies from cj = 
0.466 at T = 875 K to £, = 0.415 at T = 1250 K. An important point 
is that these values lead to similar values of the isothermal compres- 
sibility calculated using either the Carnahan Starling expression, i.e. 
pk B Tx T = (1 - 0 4 /(2£(4 - 0 + (1 - 0% or the extrapolation of 
S(q) at q = 0, i.e. pk B Tx T = S(0), S(q) being obtained from LDA 
calculations. Table I displays the close correspondence between the 



two sets of calculations, which is an indication of the thermodynamic 
consistency of our approach to determine the excess entropy. 

Using this simple analytical expression of the excess entropy, we 
plot in Fig. 7 the relationship between the excess entropy and the 
dimensionless form of the diffusion coefficient, D*. Clearly the 
Dzugutov scaling law becomes valid. We plot also the relationship 
using Eq. (10) with (d^/dT) v = 0. The results clearly indicate a 
necessity of introducing the temperature dependence of £, or a, when 
the hard sphere model is employed. 

In summary, we have computed structural and dynamic prop- 
erties of liquid aluminum by ab initio molecular dynamics using 
two different exchange and correlation potentials, LDA and GGA. 
We determine the temperature dependence of the self- diffusion coef- 
ficient as well the viscosity using direct methods. Comparisons with 
experimental data favor the LDA approximation to compute 
dynamic properties of liquid aluminum. We show that the GGA 
approximation induce more important backscattering effects due 
to an enhancement of the icosahedral short range order (ISRO) 
that impact directly dynamic properties like the self- diffusion 
coefficient. 

Our findings also support the validity of using the SE relation to 
describe dynamic properties of liquid Al on the basis of ab initio 
molecular dynamics calculations. 

Finally, we show that Dzugutov's scaling law, which relates the 
scaled diffusivity to the excess entropy, is obeyed provided the excess 
entropy as computed accurately. Our results indicate that the simple 
two -body S 2 approximation is not sufficient and that the excess 
entropy determined by the Carnahan Starling equation, including 
temperature evolution of the packing fraction, is adequate. Note that 
the excess entropy based on the Carnahan Starling equation is given 
by an analytical expression which can be easily extended to binary 
mixtures. Such a work is under consideration. 

Methods 

The AIMD simulations were performed within the framework of plane-wave-based 
density- functional theory as implemented in the Vienna ab initio simulation pack- 
age 48 . A plane- wave basis set with an energy cutoff of 241 eV was used. The electron- 
ion interaction was described by the projected augmented-wave method 49 , and both 
the local-density approximation 20 and the generalized gradient approximation 
(GGA) in the Perdew-Burke-Ernzerhof (PBE) form 21,22 for the exchange -correlation 
energy were investigated. All the dynamical simulations were carried out in the NVT 
ensemble by means of a Nose thermostat to control temperature. Newton's equations 
of motion were integrated using Verlet's algorithm in the velocity form with a time 
step of 1 fs. A cubic unit cell with periodic boundary conditions containing 256 atoms 
was used in the simulation. Only the T-point sampling was considered to sample the 
supercell Brillouin zone. The liquid samples were first prepared at 1500 K to reach 
thermal equilibrium (well above the melting temperature of 933 K) followed by 
cooling to 1250, 1 125, and 1000 K successively with a rate of 10 13 K/s. We also explore 
the undercooling state by performing runs at 875 K (LDA) and 850 K (GGA). At each 
temperature, the volume V of cell was chosen to reproduce the experimental densi- 
ties 8 and after an equilibration for 3 ps, the run was continued for 80 ps. The 
remainder of the simulation was used to extract the physical properties. 2000 con- 
figurations were used to produce averaged quantities such as the static structure factor 
and the pair-correlation function. Among these configurations, we have selected ten 
independent ones, regularly spaced in time (each separated by a time interval close to 
8 ps), to extract their inherent structures 50 . To this end, the steepest-descent energy- 
minimization procedure with the conjugate gradient method is imposed on each of 
these configurations. This method allows us to uncouple the vibrational motion from 
the underlying structural properties, since atoms are brought to a local minimum in 
the potential- energy surface. 
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